Kibble-Zurek mechanism in a quenched ferromagnetic Bose-Einstein condensate 
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The spin vortices are shown to be created through the Kibble-Zurek (KZ) mechanism in a quantum 
phase transition of a spin-1 ferromagnetic Bose-Einstein condensate, when the applied magnetic field 
is quenched below a critical value. It is shown that the magnetic correlation functions have finite 
correlation lengths, and magnetizations at widely separated positions grow in random directions, 
resulting in spin vortices. We numerically confirm the scaling law that the winding number of spin 
vortices is proportional to the square root of the length of the closed path, and for slow quench, 
proportional to Tq 1 ^ 6 with tq being the quench time. The relation between the spin conservation 
and the KZ mechanism is discussed. 
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I. INTRODUCTION 

Spontaneous symmetry breaking in a phase transition 
produces local domains of an order parameter. If do- 
mains are separated by such a long distance that they 
cannot exchange information, local domains grow ini- 
tially with random phases and eventually give rise to 
topological defects when they overlap. This scenario 
of topological-defect formation in continuous-symmetry 
breaking is known as the Kibble-Zurek (KZ) mecha- 
nism p], Q , which originally predicted the cosmic-string 
and monopole formation in the early Universe [1], and 
has since been applied to a wide variety of systems. Ex- 
perimentally, the KZ mechanism has been examined in 
liquid crystals [3|, 0], superfluid 4 He J5[ and 3 He [p. H, 
an optical Kerr medium [Hi, Josephson junctions [j|jici|], 
and superconducting films 

Recently, spontaneous magnetization in a spinor Bose- 
Einstein condensate (BEC) has attracted much interest 
as a new system for studying the KZ mechanism [l2|, EH 
EH EH- In the experiment performed by the Berkeley 
group [l2[, a BEC of F = 1 87 Rb atoms are prepared in 
the m = state, where F is the hyperfine spin and m 
is its projection on the direction of the magnetic field. 
By quench of the magnetic field, say in the z direction, 
magnetization appears in the x-y plane. Since the spinor 
Hamiltonian is axisymmetric with respect to the z axis, 
the magnetization in the x-y direction breaks the U(l) 
symmetry in the spin space. Thus, local domain forma- 
tion is expected to lead to topological defects — spin 
vortices — through the KZ mechanism. 

However, the origin of the spin vortices observed after 
the quench in the Berkeley experiment [12] cannot be at- 
tributed to the KZ mechanism. In fact, in Ref. [12], the 
spin correlation extends over the entire system (at least 
in the x direction) and the domains are not independent 
with each other. We have shown in Ref. [l3[ that the ori- 
gin of the observed spin vortices is initial spin correlation 
due to the residual m = ±1 atoms, which forms domain 
structure followed by spin vortex creation [16]. In order 



to realize the KZ mechanism in this system, i.e., in order 
to ensure that the magnetic domains grow independently, 
the size of the system must be much larger than the spin 
correlation length and the long-range correlation in the 
initial spin state must be absent. The aim of the present 
paper is to show that under these conditions spin vortices 
are generated through the KZ mechanism. 

In the present paper we will consider ID-ring and 2D- 
disk geometries. We will show that in the ID ring the 
average spin winding number after the quench is propor- 
tional to the square root of the system size, which is in 
agreement with the KZ prediction @|. In 2D the winding 
number along a path with radius R is also proportional 
to R 1 ' 2 as long as R is much larger than the vortex spac- 
ing, while it is proportional to R for small R. When 
the magnetic field is quenched slowly, the winding num- 
ber is shown to be proportional to Tq 1 ^ 6 with tq being 
the quench time. This power law can be understood by 
Zurek's simple discussion [2[. 

The spinor BEC is different from the other systems in 
which the KZ mechanism has been observed, in that the 
total spin is conserved when the quadratic Zeeman en- 
ergy q is negligible. This fact is seemingly incompatible 
with the KZ postulate, since the magnetic domains must 
be correlated with each other so that the total magnetiza- 
tion vanishes. We will show that for q = small magnetic 
domains are aligned to cancel out the local spin averaged 
over the correlation length, and that they are indepen- 
dent with each other over a greater length scale; the spin 
conservation is thus compatible with the KZ mechanism. 

The present paper is organized as follows. Section [III 
analyzes spontaneous magnetization of a spin-1 BEC and 
the resultant magnetic correlation functions using the 
Bogoliubov approximation. Section HU performs numeri- 
cal simulations of the dynamics of quenched BECs in ID 
and 2D, and shows that the KZ mechanism does emerge 
in the present system. Section HV1 provides conclusions. 
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II. BOGOLIUBOV ANALYSIS OF A 
QUENCHED FERROMAGNETIC 
BOSE-EINSTEIN CONDENSATE 

A. Hamiltonian for the spin-1 atoms 

We consider spin-1 bosonic atoms with mass M con- 
fined in a potential Vt rap (r). The noninteracting part of 
the Hamiltonian is given by 



H = J dr tin 



(r) -^ v2 + ^tra P (r-) TPm(r), 

(1) 

where ^ m (r) annihilates an atom in magnetic sublevel m 
of spin at a position r. 

The interaction between atoms with s- wave scattering 
is described by 



dr : 



cop 2 (r)+ Cl F' 2 (r) 



where the symbol :: denotes the normal order and 
1 

m= — 1 

F{r) = Y, Tpl(r)fmm4m>(r), 



(2) 

(3) 
(4) 



with / = (fx,fyifz) being the spin-1 matrices. The in- 
teraction coefficients in Eq. ([2]) are given by 



co 



Cl 



4:7rh ao + 2a2 
~M 3 ' 
Anh 2 ci2 — ao 
~M 3 ' 



(5a) 
(5b) 



where as is the s-wave scattering lengths for two colliding 
atoms with total spin S. 

When magnetic field B is applied, the linear Zeeman 
effect rotates the spin around the direction of B at the 
Larmor frequency. Since Ho and H{ nt are spin-rotation 
invariant and we assume the uniform magnetic field, the 
linear Zeeman term has only a trivial effect on spin dy- 
namics — uniform rotation of spins about B — which is 
therefore ignored. The quadratic Zeeman effects for an 
F = 1 87 Rb atom is described by 



H„ 



j4_ 

4£ h f 



m m / 



(6) 



where /ie is the Bohr magneton and E^ > is the hy- 
perfine splitting energy between F = 1 and F = 2. The 
total Hamiltonian is given by the sum of Eqs. (pQ), ([2]), 
and (O, 



B. Time evolution in the Bogoliubov 
approximation 

In the initial state, all atoms are prepared in the m = 
state. We study the spin dynamics of the system using 
the Bogoliubov approximation with respect to this initial 
state. For simplicity, we assume Vtrap = in this section. 

In the Bogoliubov approximation, the BEC part in the 
field operator is replaced by a c- number. In the present 
case, we write the m = component of the field operator 
as 



-icont/h 



(8) 



where n is the atomic density. We expand ^±i(r) as 



i>±i(r) = e 



-icont/h 



J2^=e ikr a ±1 , k , (9) 



where V is the volume of the system and a±i 5 & is the 
annihilation operator of an atom in the m = ±1 state 
with wave vector k. Keeping only up to the second order 
of S^o(r) and ^±i(r) in the Hamiltonian, we obtain the 
Heisenberg equation of motion for d±\^ as 



1 dt 



(e k +q + cin)a±i )fc (t) + ana^ _ k (i), 

(10) 

where e k = h 2 k 2 /(2M) and q = [i\B 2 / \4E hf ) . The mag- 
netic field is assumed to be applied in the z direction. The 
solution of Eq. (fTD|) is obtained as 



/.x ( E kt .e k + q + cin . E k t\ A 
fl±i,fcW = ( cos ~^ 1 ^ sin— I a±i,fc(0) 

f.an . E k t\ Af 



where 



Ek = V( £ k + q)(£k J cq J c2cin). 



(12) 



When E k is imaginary, the corresponding modes are 
dynamically unstable and grow exponentially. Since c\ < 
and q > for F = 1 87 Rb atoms, the exponential 
growth occurs for 



q < 2\a\n q c . 



(13) 



This critical value of q agrees with the phase boundary 
between the polar phase and the broken-axisymmetry 
phase [13, When q < q c /2, the wave number of 

the most unstable mode is 



, /2M fq c \ 



(14) 



H — Hq + H q + H n 



(7) and when q c /2 < q < q c , k mu = 0. 
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C. Fast quench 

We consider the situation in which q is much larger 
than the other characteristic energies for t < 0, and q is 
suddenly quenched below q c at t = 0. During t < 0, the 
time evolution in Eq. JTTJ) is a±i } fc(t) ~ e^^a^i^O), 
and the m = ±1 state remains in the vacuum state. For 
t > 0, we obtain the time evolution of the density of the 
m = ±1 component as 



sin — — 

E k ft 



1\E k \t/K 



(15) 



where the expectation value is taken with respect to the 
vacuum state of the m = ±1 component. In the sec- 
ond line of Eq. (fT5|). we have kept the unstable modes 
alone with k < k c = y/2M(q c — q) /ft by assuming that 
\Ek\t/h ^> 1. This result indicates that the m = ±1 
components grow exponentially after the quench. 

Since the operator in Eq. ([4]) is replaced by y/n in 
the Bogoliubov approximation, the magnetization oper- 



ator F+ = F_ 



- Ft _ 



iF y has the form, 



F+(r) = V2n ^{(r) + ^_i(r) 



(16) 



Using Eq. ([TTjh the time evolution of the correlation func- 
tion is calculated to be 



F+(r,t)F_(r',t) 



2n 



.Sk + q . £W 
cos — — h z — - — sm — — 



ih- (r—r) 



n ^ 
2F ^ 



2|£ fc |t/7i+ife.(r-r') 



q-Sk 



(17a) 
(17b) 



where in the second line we have kept the unstable modes 
alone. 

From the exponential factor in Eq. (|17b|) , we see that 
the sum is contributed mostly from k around the mode 
with maximum \Ek\. The denominator in the summand 
of Eq. (|17b|) is much smoother than the exponential fac- 
tor if q is not close to q c , and then we approximate Ek 
with e mu = ^ 2 /Cmu/ (2M) in the denominator. We expand 
2\Ek\t/h around k mu in the exponent as 



2\E k \t 

n 



t 



i 



eAk 2 



— S 4 Afe 4 
256 



0(Ak 6 ), 



(18) 

where Ak = k — k mu . It is clear that r sets the time 
scale for the exponential growth. The magnetization is 
observed when it sufficiently grows, i.e., t ~ r. Replacing 
the summation with the Gaussian integral in Eq. (|17b)h 



we find that £ represents the correlation length. For q < 
g c /2, k mu is given by Eq. (jHD, and 



r = 



<7c - 2g 



For (7 c /2 < q < q c , k } 



mu — and 



V M g(g c - g) " 



(19) 
(20) 

(21) 
(22) 



At (7 = g c /2, Eqs. ([20]) and ([22]) vanish, and the Afc 4 term 
in Eq. ([T8]) becomes important, with 



(— 

\2M 2 q 2 



1/4 



(23) 



We first consider a ID system with the periodic bound- 
ary condition, i.e., the ID ring geometry. We assume that 
the radius of the ring R is much larger than the domain 
size, and the curvature of the ring does not affect the 
dynamics. 

For q < q c /2, the magnetic correlation function is cal- 
culated to be 



F + (6,t)F_(6 f ,t) 



2n 



cos[k mu R(0 - 0')] 



xe 



t/r-rR 2 {e-e') 2 /{te 2 ) 



(24) 



where r and £ are given by Eqs. ([19]) and ([20]) . and 6 and 
0' are azimuthal angles. For q c /2 < q < q c: we obtain 

F + (e,t)F.(6',t)) = ±JI-*-e*'"*W % /<#) 
I Zt; y 7Tt q c — q 

(25) 

with Eqs. ((2TJ) and ([22]). At q = q c /2, the correlation 
function reads 



F+(0,t)F-(0',t) 
n q c /r\V4 



2ttE q c 



(t 



,t/r 



r| 4j oF2 U'4' ^ 



i/\4 



TR 2 (d-6')i{3\ /5 3 tR 4 (6-9') 

t — ^ — r UJ° 2 U'2' — & 



where T is the Gamma function and 

F(a)T(b) 



,(26) 



F 2 (a,b,z) = V — 
„ 1 [a 



j=0 



(a + j)T(b + j) j\ 



(27) 



is the generalized hypergeometric function. Equa- 
tion ([26]) is shown in Fig. [2](a), where H gives a char- 
acteristic width of the correlation function. 
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Next, we consider the 2D geometry. For q c /2 < q < 
g c , and then k mu = 0, the integral can be performed 
analytically, giving 



where 



UT qc c t/r-r\r-r'\ 2 /(te) 

2ir^ 2 t q c - q 



(28) 

where r and £ are given in Eqs. ([2T]) and ([22]) . For other 
g, we can perform only the angular integral as 



F + (r,t)F_(r\t) 



n q c 



47r q c - q - s n 



x / kJ (k\r-r'\)e 2 \ Ek \ t/h dk, 
Jo 

(29) 

where Jo is the Bessel function. If the exponential factor 
is much sharper than the Bessel function around /c mu , 
the correlation function ([29]) is approximated to be oc 
Mk mu \r-r f \) BE!- 

As shown above, the correlation function ()17b|) has a 
finite correlation length, and the magnetization at posi- 
tions widely separated from each other grow with inde- 
pendent directions in the x-y plane. Thus, the growth 
of the magnetic domains is expected to leave topological 
defects through the KZ mechanism. 



D. Slow quench 

In the previous sections, we have assumed that the 
magnetic field is suddenly quenched to the desired value 
at t = and q is held constant for t > 0. We assume here 
that for t > the magnetic field is gradually quenched 

as 



q(t) =q c (l- — 



The magnetic correlation can be estimated to be 



F+(r,t)F_(r',t) 

2\E k (t)\t 



oc / dkexp 



dt + ik-(r- r') 



■ (31) 



Since we are interested in the vicinity of the critical 
point where correlation starts to grow, we expand \Ek(t)\ 
around k mu = and keep the terms up to the order of 
k 2 . For the ID ring, we obtain 

(F+(e,t)F-(e',t)) oc e /(')-* 2 (*-o 2 /^ ; (32 ) 

and for the 2D geometry, 



[F+(r,t)F_(r\t)) 



J(t)-\r-r'\ 2 /e Q 



(33) 



/(*) 



TQq c 
2h 



tan 



t 



1 - 



4ft 
M 



y/t(T Q -t) 



t 

.-,1/2 



2t 



For t <C tq, f(t) can be expanded as 



/(*) 



TQq c 
2% 



s t 3 / 2 

Q 3/2 



O 



t h ' 2 

5/2 



(34) 
(35) 

(36) 



and from f(t) ~ 1, the time scale for magnetization is 
given by 



Substitution of into Eq. ([35|) yields 



\M 3 <j c 



1/6 



1/3 



(37) 



(38) 



The same power law is obtained in Ref. [14j . 

It is interesting to note that the results (|37|) and 
(f38|) are easily obtained also by the simple discussion by 
Zurek [2]. Since q(t) depends on time, r and £ given in 
Eqs. (f2Tj) and ([22]) are time dependent, and hence they 
are regarded as the growth time and correlation length 
at each instant of time. The local magnetization is de- 
veloped after a time £q has elapsed such that 



(30) Using 



r(t) 



2q c y/t(r Q -t) 2q c V~t' 
we obtain tQ in Eq. (|37|) . Substituting this into 



e(t) 



ft 2 TQ 

Mq c «(t q - i) Mq c t 



fi, 2 r Q - 2i 



(39) 



(40) 



(41) 



yields Eq. 



III. NUMERICAL SIMULATIONS AND THE 
KIBBLE-ZUREK MECHANISM 



A. Gross-Pitaevskii equation with quantum 
fluctuations 



The multicomponent Gross-Pitaevskii (GP) equation 
is obtained by replacing the field operators with the 
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macroscopic wave function ip m in the Heisenberg equa- 
tion of motion: 



ih 



±1 



dt 



2M 



V + Ftrap + q + Cop) 1p±i 



+ci ( -^=F t i/j ± F z i/j ±1 



ih 



dt 



2M 



V + Ftrap + C p tpO 



C\ 

'V2 



(42a) 



(42b) 



where p and F are defined using ijj m instead of ^ m in 
Eqs. (j3]) and The wave function is normalized as 



J dr 



(43) 



with N being the number of atoms in the condensate. 

Suppose that all atoms are initially in the m = state. 
It follows then from Eq. (|42ap that ijj±i will remain zero 
in the subsequent time evolution. This is because quan- 
tum fluctuations in the transverse magnetization that 
trigger the growth of magnetization are neglected in the 
mean-field approximation. We therefore introduce an ap- 
propriate initial noise in ip±i so that the mean- field ap- 
proximation reproduces the quantum evolution. 

Let us write the initial state as 



k vv 



(44) 



where a±\^ are c- numbers. We assume that the c- 
number amplitudes a±i are stochastic variables 
whose average values vanish, 



(a±i,fe(0))avg = 0, 



(45) 



where by (• • • ) avg we denote the statistical average over 
an appropriate probability distribution. The linear ap- 
proximation of the GP equation with respect to a±\^ 
gives the same time evolution as Eq. (fTT]L in which the 
operators are replaced by the c- numbers. We thus obtain 



0.8 
0.6 
i^0.4 
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FIG. 1: (Color online) (a) Time evolution of the auto cor- 
relation function given in Eq. (|49|) for the ID ring geometry, 
(b) Magnitude of the normalized magnetization \F+\/p (solid 
curve, left axis) and direction of the magnetization argF + 
(dashed curve, right axis) at t — 70 ms for q = and (c) for 
q = q c /2. The radius of the ring is R = 50 /xm, the atomic 
density is n — 2.8 x 10 14 cm -3 , and the number of atoms is 
iV = 10 6 . 



In the following, we will perform numerical simulation 
of spontaneous magnetization using the GP equation and 
show that the ensuing dynamics exhibits defect formation 
similar to the KZ mechanism. As the initial state of the 
m = ±1 wave functions, we use Eq. (|44j) with 



a±i,fc(0) = OL; 



rnd 



rnd i 



(48) 



where a rn d and /3 rn d are rand om v ariables following the 
normal distribution p{x) = ^2/tt exp(— 2x 2 ). Equation 
(USD then satisfies Eqs. (@5| and (|47|) . 



F+(r,t)F_(r',t) 



2n ^ 



cos — — + i—- — sm — — 

n Ek ft 



[ e -ifc-(r-r')| Olifc(0) |S 



+e 



Ik-(r-r') 



a-i,-fc(0)| 2 . 



(46) 



Comparing Eq. (|46|) with Eq. (|17a|h we find that they 
have the same form if the variance of a±i,fe(0) satisfies 



|a±i,fc(0)| 2 )av g = \ 



(47) 



for all k. 



B. ID ring geometry 

Let us first investigate the ID ring system. Experi- 
mentally this geometry can be realized, e.g. , by an opti- 
cal trap using a Laguerre-Gaussian beam [19[. We reduce 
the GP equation ([42]) to ID by assuming that the wave 
function ipm depends only on the azimuthal angle 0. The 
average density of atoms is assumed to be n = 2.8 x 10 14 
cm -3 . When the radius of the ring is R = 50 pm and 
the radius of the small circle is 2 /im, the total number 
of atoms is TV ~ 10 6 . 

Figure [1] illustrates a single run of time evolution for 
an initial state given by Eqs. and (|38|). Figure [U 
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(a) shows time evolution of the auto correlation function 
defined by 



F(t) 



I 



RdO 



\F + iP,t)? 
P 2 (0,t) 



(49) 



For both q — and q = q c /2, the transverse magneti- 
zation grows exponentially with a time constant ~ r = 
h/q c 8 ms. Snapshots of the transverse magnetization 
at t = 70 ms are shown in Figs. [T] (b) and[T] (c) for q = 
and q = q c /2, respectively. We define the spin winding 
number as 



w 



1 f 2lT 



Rd6 



1 



2i\F + \ 



dF + 
~d0~ 



dF_ 
89 



(50) 



which represents the number of rotation of the spin vector 
in the x-y plane along the circumference, and of course 
w is an integer. The spin winding numbers are w = 7 in 
Fig. □ (b) and w = -1 in Fig. [Q (c). 

Figure [2] (a) shows the ensemble average of the nor- 
malized correlation function, 



(F corr (58)) 



avg 



JdOF + (O)F_(9 + S0)\ 
J d0p(6)p(0 + 50) / 



(51) 



avg 



at t = 70 ms. For q = q c /2, the correlation function has 
the characteristic width of ~ 5 in Eq. ([23]) , indicating 
that the ring is filled with magnetic domains with an av- 
erage size of ~ S. According to the KZ theory, the mag- 
netic domains with random directions give rise to the spin 
winding, which is estimated to be w ~ (R/S) 1 / 2 . This R 
dependence of w is clearly seen in Fig. [51(b). The ensem- 
ble average of the winding number, (w) avg , vanishes due 
to the random nature of the initial noise, and the square 
root of its variance, (w 2 ) a vg, should be regarded as a typ- 
ical winding number. The variance is expected to obey 
the x 2 distribution with 1000 degrees of freedom, and 
hence we show the 95% confidence interval to estimate 
the statistical errors in Fig. [2j As shown in the inset of 
Fig. [2] (b), the typical winding number changes in time, 
since the ferromagnetic energy is converted to the kinetic 
energy and the system exhibits complicated dynamics. 

The situation is different for q = 0, in which the cor- 
relation function oscillates with a Gaussian envelope as 
shown in Fig. [2] (a). This form of the correlation function 
gives us the answer to the question as to how the KZ 
mechanism manifests itself in spin conserving systems. 
The finite correlation length for q = indicates that the 
spin is conserved not only globally but also locally, that 
is, the locally integrated spin over the correlation length 



\Sr\<£ 



F(r + Sr)dSr, 



(52) 



is held to be zero for any r. This local spin conservation 
is due to formation of staggered domain or helical spin 
structures whose periodic length is much smaller than 



(a) 0.2 



f 0.1 



05 



-0.1 k 

-71/1 
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q=0 | A 











58 



n/1 



(b) 1000 
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FIG. 2: (Color online) (a) Numerically obtained correlation 
function given in Eq. (|51|) at t = 70 ms (solid curves), and 
theoretical fits (dashed curves) from Eqs. ([24]) and ([26]) . Other 
parameters are the same as those in Fig. [1] (b) R dependence 
of the variance of the spin winding number, where the number 
of atoms is related to R as N = 10 6 x R [fim] /50. The dashed 
lines are semi-log fits to the numerical data. The inset shows 
the time dependence of {w 2 ) avg for R = 50 fim. The data 
in (a) and (b) are averages over 1000 runs of simulations for 
different initial states produced by random numbers. The 
error bars in (b) represent the 95% confidence interval of the 
X 2 distribution. 



£. Thus, the neighboring domains tend to have opposite 
magnetizations to cancel out the spin locally, and the 
domains far from each other grow independently; the spin 
conservation and the KZ mechanism are thus compatible. 

The oscillation in the correlation function originates 
from the fact that the most unstable modes have nonzero 
wave numbers ±fc mu . Each correlated region of size ~ 
£ = [8h 2 /(Mqc)} 1 / 2 contains spin waves of e lkmuRe and 
e~ lkmuRe . If there is an imbalance between these modes, 
the winding number monotonically increases or decreases 
in each region of ~ £. This is the reason why (w 2 ) dcvg is 
larger for q = than for q = q c /2 in Fig.[2](b). It follows 
from this consideration that for k mu ^ ^> 1 the winding 
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100 




10 



l-2q/q c 



FIG. 3: (Color online) Dependence of the variance of the 
spin winding number on q. Except for q, the parameters are 
the same as those in Fig. [T] The dashed line is proportional 
to (1 — 2q/q c ) 3 / 2 . The plots show the averages over 1000 runs 
of simulations for different initial states produced by random 
numbers. The error bars represent the 95% confidence inter- 
val of the x 2 distribution. 



number is proportional to 



w ~ k mu £ x — = k n 



R£ oc ( 1 - ^ 



3/4 



(53) 



where Eqs. (p~4|) and (|20|) are used. Figure [3] shows the 
averaged variance of the winding number versus 1 — 2q/q c . 
For small g, (w 2 ) avg is proportional to (1 — 2q/q c ) 3 / 2 , in 
agreement with Eq. ([53]) . When q is close to q c /2, the 
spin winding within the correlated region, fc mu £, becomes 
small, and then the winding number reduces to the value 
shown in Fig. [21(b), i.e., {w 2 ) avg — 4. 

We next discuss the results of simulations of slow 
quench as in Eq. (|30|) . Since the winding number for 
the slow quench is small compared with the fast quench, 
we take a large ring of R = 400 jam. Figure 0] shows 
the variance of the winding number as a function of the 
quench time. We can clearly see that (w 2 ) SiVg has a power 

— 1/3 

law of Tq within the statistical error, which is in agree- 



ment with ^q 1 ~ Tq , with £q being given in Eq. 
Thus, the present system follow the quench- time scaling 
of Zurek [2j. We note that, in the slow quench, the wind- 
ing number converges to an almost constant value for 
varying quench time tq, as shown in the inset of Fig. [H 
This is because little excess energy other than for exciting 
spin vortices is available for the slow quench. 



C. 2D disk geometry 

When the confinement in the z direction is tight, the 
system is effectively 2D. For simplicity, we ignore the den- 
sity dependence in the z direction, and assume that the 
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FIG. 4: (Color online) Dependence of the variance of the spin 
winding number at t = 400 ms on the quench time tq , where q 
is varied as in Eq. ([30]) . The radius of the ring is R = 400 /im, 



the atomic density is n — 
of atoms is N = 8 x 10 6 



tO T, 



-1/3 



2.8 x 10 cm , and the number 
The dashed line is proportional 
The inset shows time evolution of {w 2 ) avg . The 
data are averages over 1000 runs of simulations for different 
initial states produced by random numbers. The error bars 
represent the 95% confidence interval of the \ 2 distribution. 



2D GP equation has the same form as Eq. (|^2]) . We as- 
sume that the wave function vanishes at the wall located 
at (x 2 + y 2 ) 1 / 2 = R w = 100 /im, and that the potential 
is flat inside of the wall. Then the density n = 2.8 x 10 14 
cm -3 is almost constant except within the healing length 
{3/[87rn(a + 2a 2 )]} 1/2 ~ 0.16 /im near the wall. When 
the thickness in the z direction is ^ 1 /im, the number of 
atoms is N ~ 10 T . Such a system will be realized using 
an optical sheet and a hollow laser beam. 

The initial state of tpo is a stationary solution of the GP 
equation, and the initial state of i/j±i is given by Eq. (|44j) 
with random variables flU]) . Figure [5] (a) shows time 
evolution of the autocorrelation function of the transverse 
magnetization, 



F(t) 



dr 



\F+(r,t)f 
P 2 (r,t) 



(54) 



which grows exponentially with the same time constant 
as that in Fig. [TJ and saturates for t > 100 ms. 

Snapshots of |F + (r)| and argF + (r) at t = 100 ms are 
shown in Figs. [5] (b) and [5] (c). We see that |F+(r)| 
at t ^ 100 ms contains many holes, around which the 
spin direction rotates by 2tt. Since this topological spin 
structure consists of singly-quantized vortices in the m = 
±1 states filled by atoms in the m = state, it is called 
the "polar-core vortex." We can estimate the spin healing 
length £ s by equating the kinetic energy ft 2 /(2M^ 2 ) with 
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t = 100 ms 
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(c) q = q / 2 t = 50 



FIG. 5: (Color) (a) Time evolution of the autocorrelation function given in Eq. ([54]) for the 2D disk geometry. The radius of 
the disk is R w — 100 /mi, the atomic density is n — 2.8 x 10 14 cm -3 , and the number of atoms is N — 10 7 . (b) Profiles of the 
magnetization (upper) and its direction argF + (lower) for q = and (c) for q = q c /2. The size of each panel is 200 /im 
x 200 /mi. 



the energy of magnetization \q — q c \, giving 



^2M\q-q c \ 

This length scale is £ s c± 1.7 /im for q — and £ s c± 2.4 
/im for g = <7c/2, which are in good agreement with the 
sizes of the vortex cores in Figs. [5] (b) and [5] (c). 



In 2D, the correlation function is defined by 

iw (K u / J drF + (r)F_(r + 5r) \ 
( F ^(Sr))^ = { Jdrp{r)p{r + Sr) , (56) 

which are shown in Figs. [6] (a) and[6](b). We find that as 
in the ID case the most unstable wave length is reflected 
in the shape of the spin correlation function ([56]) , and the 
characteristics of these correlation functions in the ra- 
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FIG. 6: (Color) (a) Spin correlation function denned in 
Eq. (|56l) at t = 100 ms for q = and (b) for g g c /2. 
(c) The variance of the winding number along the circumfer- 
ence of the circle of radius R. The dashed lines and dotted 
lines are proportional to R and R 2 , respectively. In (a)-(c) 
the parameters are the same as those in Fig. [5] and the data 
are averages over 1000 runs of simulations for different initial 
states produced by random numbers. The error bars in (c) 
represent the 95% confidence interval of the % 2 distribution. 



dial direction are similar to those in ID shown in Fig. [2j 
For q = 0, the mean distance between spin vortices in 
Fig. [5] (b) is not determined by the correlation length 
(the whole width of the concentric pattern in Fig. [6] (a)) 
but by ~ fc~J, i.e., the width of the concentric rings in 
Fig. [6] (a). On the other hand, for q = q c /2, the density 
of spin vortices is determined by the correlation length, 
i.e., the size of the blue circle ~ 30 jam in Fig.[6](b). The 
staggered concentric correlation for q = suggests that 
the spin is conserved locally within the region of the cor- 
relation length, and domains at a distance larger than the 
correlation length grow independently, while conserving 
the total spin. 

The spin winding number for 2D is defined as 



w(R) = 



2tt 



1 



C(R) 



2i\F+ 



(F_VF+ -F+VF-) -dr, 



(57) 

where C(R) is a circle with radius R < R w located at the 
center of the system. Figure [6] (c) shows the R depen- 



dence of the ensemble average of w 2 (R), where the radius 
of the system is fixed to R w = 100 /im and the data are 
taken at t = 100 ms. It should be noted that (w 2 (R)) diVg 
is proportional to R for large i?, as expected from the 
KZ theory [2|, while it is proportional to R 2 for small R. 
This R 2 dependence is due to the fact that the probabil- 
ity P for a spin vortex to be in the circle is proportional 
to ttR 2 . The variance of the winding number is therefore 
0(l-P) + l 2 P/2 + (-l) 2 P/2 ex R 2 , if the probability that 
two or more vortices enter the circle is negligible. This 
condition is met when the density of spin vortices times 
ttR 2 is much smaller than unity, and hence the radius 
R at which the crossover from (w 2 (R)) 8CVg oc R to oc R 2 
occurs is larger for q = q c /2 than for q = 0. As in ID, 
nonzero k mu enhances the winding of magnetization, and 
the winding number is larger for q = than for q = q c /2. 

Figures [5] (b) and [51(c) obviously show that the density 
of spin vortices is uniform when the size of the system is 
large enough. The number of spin vortices in a radius R 
is therefore proportional to R 2 . If the topological charge 
of each spin vortex, +1 or —1, was chosen at random, 
the net winding number along the circle of radius R, i.e., 
the difference between the numbers of +1 and —1 vor- 
tices would be proportional to R. However, from Fig. [6] 
(c), the winding number is proportional to R}l 2 for large 
P, consistent with the KZ mechanism. The topological 
charge of each spin vortex is thus not at random but 
anticorrelated to each other to reduce the net winding 
number. 

Figure [71 shows the result of the slow quench for 2D, 
where q(t) is given by Eq. (|30|) . The winding number fol- 
lows the scaling law, (w 2 ) avg oc Tq 1 ^ 3 , as predicted from 
Eq. (|38]h indicating that Zurek's discussion is applicable 
also to 2D. In order to obtain this scaling law, we must 
specify the time at which the winding number is taken, 
since the spin winding number decays in time, as shown 
in the inset of Fig. From the scaling law in Eq. (|3T|h 
we specify the time to take the winding number as 



t 



1/3 



(I) 



2/3 



const., 



(58) 



which is indicated by the arrows in the inset of Fig. [3 



IV. CONCLUSIONS 

In this paper, we have studied the dynamics of a spin-1 
BEC with a ferromagnetic interaction after quench of the 
applied magnetic field in an attempt to investigate spon- 
taneous defect formation in the spinor BEC. We have 
analyzed the magnetization triggered by quantum fluc- 
tuations using the Bogoliubov approximation, and per- 
formed numerical simulations of the GP equation with 
initial conditions that simulate quantum fluctuations. 

We have shown that the correlation functions of the 
magnetization have finite correlation lengths (Figs. [2j [6] 
(a), and[6](b)), and therefore magnetic domains far from 
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FIG. 7: (Color online) (a) Variance of the spin winding num- 
ber versus the quench time tq for the 2D disk geometry, where 
q is varied as in Eq. ([30]) . The inset shows time evolution of 
{w 2 ) avg . The plots are taken at the times when t/r^ 3 = 
constant is satisfied, which are shown by the arrows in the 

2 /3 

inset. The dashed line is proportional to Tq . The radius 
of the disk is R w = 400 /im and the closed path for taking 
the winding number is R = 320 /im. The atomic density is 
n — 2.8 x 10 14 cm -3 and the number of atoms is N — 1.6 x 10 8 . 
The data are averages over 1000 runs of simulations for dif- 
ferent initial states produced by random numbers. The error 
bars represent the 95% confidence interval of the % 2 distribu- 
tion. 



each other grow in random directions. We find that topo- 
logical defects — spin vortices — emerge through the KZ 
mechanism. We have confirmed that the winding num- 
ber along the closed path is proportional to the square 
root of the length of the path (Figs. [2] (b) and [6] (c)), 
indicating that the topological defects are formed from 
domains with random directions of magnetizations. 

Even when the total magnetization is conserved for 
q = 0, the winding number has the same dependence on 
the length of the path (Fig. [2] (b)). This is due to the 
fact that domains within the correlation length tend to 
be aligned in such a manner as to cancel out local mag- 
netization, and consequently the total magnetization is 
conserved. Thus, the neighboring domains have local cor- 



relation, while domains far from each other are indepen- 
dent, which makes the KZ mechanism compatible with 
the total spin conservation. The formation of the local 
correlation also creates topological defects as well as the 
KZ mechanism, and the winding number exhibits the q 
dependence as shown in Fig. [3j 

When the magnetic field is quenched in finite time tq 
as in Eq. ([30]) , the winding number has been shown to 

be proportional to Tq 1 ^ 6 (Figs. [Hand [7j). This tq depen- 
dence of the winding number can be understood from 
Zurek's simple discussion [2]: the domains are frozen at 
which the spin relaxation time becomes the same order 
of elapsed time. 

In the Berkeley experiment [l2| , the system is an elon- 
gated quasi-2D geometry, and not suitable for testing 
the KZ mechanism. The KZ mechanism should apply 
to the system in which the size of the system in the x 
direction is made much larger. In this case, the har- 
monic potential may affect the scaling law, which mer- 
its further study. Moreover in the experiment, from the 
analysis in Ref. [l3j], there are some initial atoms in the 
m = ±1 components with long-range correlation, which 
play a role of seeds for large domains and hinder the ob- 
servation of the KZ mechanism. If the residual atoms in 
the m = ±1 components is eliminated completely, mag- 
netization is triggered by quantum fluctuations as shown 
in the present paper. Another way to remove the effect 
of the residual atoms may be applying random phases to 
the m = ±1 states to erase the initial correlation. 

Note added. After our work was completed, the 
preprint by Damski and Zurek [20[ appeared, which per- 
forms ID simulations of the quench dynamics of a spin-1 
BEC. 
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